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Complex fluids exhibit structure on a wide range of length and time scales, and 
hierarchical approaches are necessary to investigate all facets of their often unusual 
properties. The study of idealized coarse-grained models at different levels of coarse- 
graining can provide insight into generic structures and basic dynamical processes at 
equilibrium and non-equilibrium. 

In the first part of this lecture, some popular coarse-grained models for membranes 
and membrane systems are reviewed. Special focus is given to bead-spring models with 
different solvent representations, and to random-interface models. Selected examples of 
simulations at the molecular and the mesoscopic level are presented, and it is shown how 
simulations of molecular coarse-grained models can bridge between different levels. 

The second part addresses simulation methods for complex fluids under shear. After 
a brief introduction into the phenomenology (in particular for liquid crystals), different 
non-equilibrium molecular dynamics (NEMD) methods are introduced and compared 
to one another. Application examples include the behavior of liquid crystal interfaces 
and lamellar surfactant phases under shear. Finally, mesoscopic simulation approaches 
for liquid crystals under shear are briefly discussed. 

1 Introduction 

The term "Complex fluid" or "soft condensed matter" - these two are often used interchangeably 
- refers to a broad class of materials, which are usually made of large organic molecules and have 
a number of common features [fl~||2]: They display structure on a nanoscopic or mesoscopic scale; 
characteristic energies are of the order of fcgT at room temperature, hence the properties are to 
a large extent dominated by entropic effects, and the materials respond strongly to weak external 
forces [3|; the characteristic response times span one or several orders of magnitude, and the The- 
ological properties of the fluids are typically non-Newtonian J2]|4]. Some examples are polymer 
melts and solutions [5 6,|7]|8), emulsions, colloids |9] [TO] QT|, amphiphilic systems lfl2l [T3l . and 
liquid crystals 01 [B] [161. 

Computer simulations of complex fluids are particularly challenging, due to the hierarchy of 
length and time scales that contributes to the material properties. With current computer resources, 
it is practically impossible to describe all aspects of such a material within one single theoretical 
model. Therefore, one commonly uses different models for different length and time scales. The 
idea is to eliminate the small-scale degrees of freedom in the large-scale models, and to incorporate 
the details of the small-scale properties in the parameters of a new, simpler, model. This is called 
coarse-graining. 

There are two aspects to coarse-graining. First, systematic coarse-graining procedures must be 
developed and applied in order to study materials quantitatively on several length and time scales. 
This is the challenge of multiscale modeling, one major growth area in materials science. It is, 
however, not the subject of the present lecture. Some recent reviews can be found in Refs. |[T7l[T8l 

Second, coarse-grained idealized models are used to study generic properties of materials on a 
given scale, and physical phenomena which are characteristic for a whole class of materials. Here, 
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microscopic details are disregarded because they are deemed irrelevant for the physical properties 
under consideration. This approach to materials science has a long tradition, going back to Ising's 
famous lattice model for magnetism. 

In the present lecture, we shall concentrate on coarse-grained models of the second type. We shall 
discuss how such models can be constructed, how they can be studied by computer simulations, and 
how such studies help to improve our understanding of equilibrium and nonequilibrium properties 
of complex materials. We will not attempt to present an overview - this would require a separate 
book - but rather focus on selected case studies. In the first part, we will consider coarse-grained 
model systems for amphiphilic membranes and their use for studies of equilibrium properties. In the 
second part, we will turn to complex fluids under shear, i.e., nonequilibrium systems, and discuss 
coarse-grained approaches for membrane systems and liquid crystals. 

2 Coarse-grained models for surfactant layers 

In this section, we shall discuss coarse-grained simulation models and methods for amphiphilic 
membranes and membrane stacks. After a brief introduction into the phenomenology and some 
general remarks on coarse-graining, we will focus on two particular types of coarse-grained mod- 
els: Particle-based bead-spring models and mesoscopic membrane models. We shall present some 
typical examples, and illustrate their use with applications from the literature. 

2.1 Introduction 

Amphiphilic membranes are a particular important class of soft material structures, because of their 
significance for biology fl9l l20l |2D . Biomembranes are omnipresent in all living beings, they 
delimit cells and create compartments, and participate in almost all biological functions. Figure Q] 
shows an artist's view on such a biomembrane. It is mainly made of two coupled layers of lipids, 
which serve as a matrix for embedded proteins and other molecules. 



This structure rests on a propensity of lipids to self-assemble into bilayers when exposed to 
water. The crucial property of lipids which is responsible for this behavior is their amphiphilic 
character, i.e., the fact that they contain both hydrophilic (water loving) and hydrophobic (water 
hating) parts. Most lipid molecules have two non polar (hydrophobic) hydrocarbon chain (tails), 




Figure 1: Artist's view of a biomembrane. Source: NIST 
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which are attached to one polar (hydrophilic) head group. Figure [2]l shows some typical examples. 
The tails vary in their length, and in the number and position of double bonds between carbon atoms. 
Whereas single bonds in a hydrocarbon chain are highly flexible, in the sense that chains can rotate 
easily around such a bond, double bonds are stiff and may enforce kinks in the chain. Chains with 
no double bonds are called saturated. The variety of polar head groups is even larger. Overviews can 
be found in |fl9l . 

By assembling into bilayers, the lipids can arrange themselves such that the head groups shield 
the tails from the water environment. Therefore, lipids in lipid-water mixtures often tend to form 
lamellar stacks of membranes separated by thin water layers (Fig. [3j- Such stacks can typically hold 
up to 20 % water. In the more dilute regime, membranes may close up to form vesicles or other 
more disordered structures. Vesicles play an important role in biological systems, since they provide 
closed compartments which can be used to store or transport substances. Single, isolated membranes 
are metastable with respect to dissociation for entropic reasons, but they have very long lifetimes. 

Depending on the particular choice of lipids and the properties of the surrounding aqueous fluid 
(pH, salt content etc.), several other structures can also be observed in lipid-water systems - spherical 
and cylindrical micelles, ordered structures with hexagonal or cubic symmetry etc. El . In this 
lecture, we shall only regard bilayer structures. 

Membranes also undergo internal phase transitions. Figure|4]shows schematically some charac- 
teristic phases of one-component membranes. The most prominent phase transition is the "main" 
transition, the transition from the liquid state into one of the gel states, where the layer thickness, 
the lipid mobility, and the conformational order of the chains, jump discontinuously as a function 
of temperature. In living organisms, most membranes are maintained in the liquid state. Neverthe- 
less, the main transition presumably has some biological relevance, since it is very close to the body 
temperature for some of the most common lipids (e.g., saturated phospholipids) lfl9l . 

In sum, several aspects of membranes must be studied in order to understand their structure and 
function, in biological as well as in artificial systems: 

• Self-assembling mechanisms 

• The inner structure and internal phase transitions, and the resulting membrane properties (per- 
meability, viscosity, stiffness) 

• The organization of membranes on a mesoscopic scale (vesicles, stacks, etc.) 

• The interaction of membranes with macromolecules. 



Figure 2: Examples of membrane lipids. Left: Phosphatidylcholine (unsaturated); Right: Sphin- 
gomyelin. 
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Figure 3 



Self-assembled bilayer structures. 



These topics relate to the physics and chemistry of amphiphiles on many length and time scales, 
and several driving forces contribute to the unique properties of membranes: On the scale of atoms 
and small molecules, one has the forces which are responsible for the self-assembly: The hydropho- 
bic effect, and the interactions between water and the hydrophilic head groups. On the molecular 
scale, one has the interplay of local chain conformations and chain packing, which is responsible 
for the main transition and determines the local structure and fluidity of the membrane. On the 
supramolecular scale, one has the forces which govern the mesoscopic structure and dynamics of 
membrane systems: The membrane elasticity, the thermal fluctuations, the hydrodynamic interac- 
tions with the surrounding fluid, the thermodynamic forces that drive phase separation in mixed 
membranes etc. 

Studying all these aspects within one single computer simulation model is clearly impossible, 
and different models have been devised to study membranes at different levels of coarse-graining. 
Among the oldest approaches of this kind are the Doniach [23] and the Pink model ll24l|25ll , gen- 
eralized Ising models, which have been used to investigate almost all aspects of membrane physics 
on the supramolecular scale ||2T1|261 . With increasing computer power, numerous other simulation 
models have become treatable in the last decades, which range from atomistic models, molecular 
coarse-grained models, up to mesoscopic models that are based on continuum theories for mem- 
branes. 

In the following, we shall focus on two important examples: Coarse-grained chain models de- 
signed to describe membranes on the molecular scale, and random interface models that describe 
membranes on the mesoscopic scale. We stress that this our overview is not complete, and there 
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Figure 4: Some prominent membrane phases. 



4 



exist other successful approaches, which are not included in our presentation. A relatively recent 
review can be found in Ref . Il27l . 

2.2 Coarse-grained molecular models 

In coarse-grained molecular models, molecules are still treated as individual particles, but their struc- 
ture is highly simplified. Only the ingredients deemed essential for the material properties are kept. 
For amphiphiles, these are the amphiphilic character (for the self-assembly), and the molecular flex- 
ibility (for the internal phase transitions). 

2.2.1 Spring-bead models 

The first coarse grained molecular models have been formulated on a lattice, relatively shortly after 
the introduction of the above-mentioned Pink model [ 24 25 1 . One particularly popular model is the 
Larson model [28 , 29 1: Water molecules (w) occupy single sites of a cubic lattice, and amphiphile 
molecules are represented by chains of "tail" (t) and "head" (h) monomers. Only particles on neigh- 
bor lattice sites interact with each another. The lattice is entirely filled with w, h, or t particles. If 
one further takes all hydrophilic particles, h and w, to be identical, the interaction energy is deter- 
mined by a single interaction parameter, which describes the relative repulsion between hydrophobic 
and hydrophilic particles. The model reproduces self-assembly and exhibits many experimentally 
observed phases, i.e., the lamellar phase, the hexagonal phase, the cubic micellar phase, and even 
the bicontinuous gyroid phase l29l [30ll . 

Lattice models can be simulated efficiently by Monte Carlo methods. However, they have the 
obvious drawback of imposing an a priori anisotropy on space, which restricts their versatility. 
For example, internal membrane phase transitions and tilted phases cannot be studied easily. In 
dynamical studies, one is restricted to using Monte Carlo dynamics, which is unrealistic. Therefore, 
most more recent investigations rely on off-lattice models. 

Off-lattice molecular approaches to study self-assembling amphiphilic systems have been ap- 
plied for roughly 15 years. Presumably the first model was introduced by B. Smit et al. ||3T1 in 1990. 
A schematic sketch is shown in Fig.|5]a). As in the Larson model, the amphiphiles are represented by 
chains made of very simple h- or f-units, which are in this case spherical beads. "Water" molecules 
are represented by free beads. Beads in a chain are joined together by harmonic springs, where a 
cutoff is sometimes imposed in order to ensure that the beads cannot move arbitrarily far apart from 
one another. Non-bonded beads interact via simple short-range pairwise potentials, e.g., a truncated 
Lennard-Jones potential. The parameters of the potentials are chosen such that ht pairs and hw 
pairs effectively repel each other. The Smit model reproduces self-assembly, micelle formation and 
membrane formation 11321 l33l l34l 1351 l36ll . 

Most molecular off-lattice models are modifications of the Smit model |27l [37 38, 39]. They 
differ from one another in the specific choice of the interaction potentials between non-bonded beads, 
in the bond potentials between adjacent beads on the chain, in the chain architecture (number of 
head and tail beads, number of tails), and sometimes contain additional potentials such as bending 
potentials. Nowadays, Smit models are often studied with dissipative particle (DPD) dynamics, and 
the interaction potentials are the typical soft DPD potentials ll37l [39l . A "water particle" is then 
taken to represent a lump of several water molecules. An introduction into the DPD method is given 
in B. Dunweg's lecture (this book). It works very well for membrane simulations as long as the time 
step is not too large |40|. 

This approach is presented in more detail in B. Smit's seminar (this book). 
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2.2.2 Implicit solvent and phantom solvent 

The primary requirement for a lipid membrane model is to ensure that the lipids maintain a bilayer 
structure. In the models discussed in the previous section [2.2.11 the bilayers are stabilized by the 
surrounding solvent particles. If the interactions between w, h, and t particles are chosen suitably, the 
solvent particles and the amphiphiles arrange themselves such that the hydrophilic and hydrophobic 
particles are separated (microphase separation). 

This "explicit solvent" approach is in some sense the most natural one and has several advan- 
tages: The bilayers are fully flexible, and they are surrounded by a fluid. The fact that the fluid 
is composed of large beads as opposed to small water molecules is slightly problematic, since the 
solvent structure may introduce correlations, e.g., in simulations with periodic boundary conditions. 
However, the correlations disappear if the membranes or their periodic images are separated by a 
large amount of solvent. 

On the other hand, explicit solvent models also have a disadvantage: The simulation time de- 
pends linearly on the total number of beads. A large amount of computing time is therefore spent on 
the uninteresting solvent particles. Therefore, effort has been spent on designing membrane models 
without explicit solvent particles. 

One early idea has been to force the amphiphiles into sheets by tethering the head groups to two- 
dimensional opposing surfaces, e.g., by a harmonic potential, or by rigid constraints BT1 l42l l43l l44l 
[45| (Fig.[5]b)). This model, which we call "sandwich model", is simple, extremely efficient, and the 
configurations are easy to analyze, since the bilayer is always well-defined. On the other hand, the 
bilayer has no flexibility, i.e., undulations and protrusions |46] are supressed and important physics 
is lost. 

A second approach is to eliminate the solvent, and represent its effect on the amphiphiles by 
appropriate effective interactions between monomers (Fig. c)). This way of modeling solvent is 
common in polymer simulations. Nevertheless, producing something as complex as membrane self- 
assembly is a non-trivial task. The first implicit solvent model was proposed only rather recently, 
in 2001, by Noguchi and Takasu ll47ll . They mimick the effect of the solvent by a multibody "hy- 
drophobic potential", which depends on the local density of hydrophobic particles. From a physical 
point of view, the introduction of multibody potentials in an implicit solvent model is reasonable 
- multibody potentials emerge automatically whenever degrees of freedom are integrated out sys- 
tematically. Nevertheless, people usually tend to favor pair potentials. Farago |48ll and Cooke et 
al. |49| have shown that it is possible to stabilize fluid bilayers with purely pairwise interactions in a 
solvent-free model. The crucial requirement seems to be that the interactions between hydrophobic 
beads have a smooth attractive part of relatively long range |49l . 

The clear benefit of implicit solvent models is their high efficiency. They are very well suited 
to study large scale problems like vesicle formation ll47l l49l . vesicle fusion [50|, vesicles subject 
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Figure 5: Schematic picture of bead-spring models for self-assembling membranes, (a) Explicit 
solvent model (b) Sandwich model (c) Solvent free model (d) Phantom solvent model. 
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to external forces BTl . phase separation on vesicles [49 1 etc. On the other hand, they also have 
limitations: They require tricky potentials; the pressure is essentially zero (since the membranes are 
surrounded by empty space); also, they cannot be used for dynamical studies where the hydrody- 
namical coupling with the solvent becomes important. 

A third approach was recently put forward by ourselves: the "phantom solvent" model |45l 
(Fig.[5Jl)). We introduce explicit solvent particles, which do not interact with one another, only with 
the amphiphiles. The amphiphiles perceive them as repulsive soft beads. In the bulk, they simply 
form an ideal gas. They have the simple physical interpretation that they probe the free volume which 
is still accessible to the solvent, once the membrane has self-assembled. Thus the self-assembly is 
to a large extent driven entropically. 

In Monte Carlo simulations, this model is just as efficient as solvent-free models, because Monte 
Carlo moves of phantom beads that are not in contact with a membrane are practically free of (com- 
putational) charge. It is robust, we do not need to tune the lipid potentials in order to achieve self- 
assembly. Pressure can be applied if necessary, and with a suitable dynamical model for the solvent 
(e.g., dissipative particle dynamics), we can also study (hydro)dynamical phenomena. Compared 
with Smit-type models, the phantom solvent model has the advantage that the solvent is structure- 
less and cannot induce artificial correlations. On the other hand, the fact that it is compressible like 
a gas, may cause problems in certain dynamical studies. 

To summarize, we have presented several possible ways to model the solvent environment of a 
self-assembled membrane, and discussed the advantages and disadvantages. The different models 
are sketched schematically in Fig. [5] 

We shall now illustrate the potential of molecular coarse-grained models with two examples. 

2.2.3 First application example: The ripple phase 

Our first example shows how coarse-grained simulations can help to understand the inner structure of 
membranes. We have introduced earlier the liquid and gel phases in membranes (Fig. [4]). However, 
one rather mysterious phase was missing in our list: The ripple phase (Pp>). It emerges as an 
intermediate state between the tilted gel phase (Lp>) and the liquid phase (L a ), and according to 
freeze-etch micrographs, AFM pictures, and X-ray diffraction studies, it is characterized by periodic 
wavy membrane undulations. More precisely, one observes two different rippled structures: Figure|6] 
shows electron density maps of these two states, as calculated from X-ray data by Sengupta et 
al 11521 . The more commonly reported phase is the asymmetric ripple phase, which is characterized 
by sawtooth profile and has a wavelength of about 150 Angstrom. The symmetric rippled phase 
tends to form upon cooling from the liquid L a phase and has a period twice as long. It is believed 
to be metastable. The molecular structure of both rippled states has remained unclear until very 
recently. 

Atomistic simulations of de Vries et al, have finally shed light on the microscopic structure of 
the asymmetric ripple phase l53l . These authors observed an asymmetric ripple in a Lecithin bilayer, 
which had a structure that was very different from the numerous structures proposed earlier. Most 
strikingly, the rippled membrane is not a continuous bilayer. The ripple is a membrane defect where 
a single layer spans the membrane, connecting the two opposing sides (Fig. [7j». 

The simulations of de Vries et al seem to clarify the structure of the asymmetric ripple phase, 
at least for the case of Lecithin. Still the question remains whether their result can be generalized 
to other lipid layers, or whether it is related to the specific properties of Lecithin head groups. This 
question can be addressed with idealized coarse-grained models. 

Rippled phases had also been observed in simulations of coarse-grained models. Kranenburg et 
al have reported the existence of a rippled state in a Smit model with soft DPD interactions B4ll55ll . 
The structure of their ripple is different from that of Fig. [7] Judging from this result, one is tempted 
to conclude that the structure of de Vries et al may not be generic. 
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Figure 6: Electron density maps of (top) an asymmetric rippled state (DMPC at 18.2 degrees) and 
(bottom) a symmetric rippled state (DPPC at 39.2 °C). Reprinted with permission from Sengupta, 
Raghunathan, Katsaras ll52l (rescaled). Copyright 2003 by the American Physical Society. 



However, we recently found ripples with a structure very similar to that of Fig. [7] in simulations 
of a simple bead-spring model by (Lenz and Schmid [56, 57|). Our lipids are represented by linear 
chains with six tail beads (size a) and one head bead (size 1.1 a), which are connected by springs 
of size 0.7cr. Tail beads attract one another with a truncated and shifted Lennard-Jones potential. 
Head beads are slightly larger and purely repulsive. The solvent environment is modelled by a fluid 
of phantom beads of size a. The system was studied with Monte Carlo simulations at constant 
pressure, using a simulation box of fluctuating size and shape in order to ensure that the pressure 
tensor is diagonal. 

Membranes were found to self-assemble spontaneously at sufficiently low temperatures, they 
exhibit a fluid L Q -phase as well as a tilted gel L^-phase. The two phases are separated by a temper- 
ature region where ripples form spontaneously. In small systems, the ripples are always asymmetric. 
In larger systems, asymmetric ripples form upon heating from the Lg/ -phase. Upon cooling from 
the L Q -phase, a second, symmetric, type of ripple forms, which has twice the period than the asym- 
metric ripple and some structural similarity, except that the bilayer remains continuous (see Fig . [8} . 




Figure 7: Asymmetric ripple in a Lecithin Bilayer, as observed in an atomistic simulation. From de 
Vries et al 11551 . Copyright 2005 National Academy of Sciences, U.S.A. 
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Figure 8: Ripples from simulations of a coarse-grained model, (a) Asymmetric ripple (b) Symmetric 
ripple. 

The sideview of this second ripple has striking similarity with the density maps of Fig.|6]b). 

The coarse-grained simulations of Lenz thus not only demonstrate that the asymmetric ripple 
structure found by de Vries is generic, and can be observed in highly idealized membrane models as 
well. They also suggest a structural model for the yet unresolved symmetric ripple. 

2.2.4 Second application example: Membrane stacks 

In our second example, we consider the fluctuations and defects in stacks of membranes. We shall 
show how simulations of coarse-grained membrane models can be used to test and verify mesoscopic 
models, and how they can bridge between coarse-graining levels. 

At the mesoscopic level, membranes are often represented by random interfaces (see also Sec- 
tion |23]l. One of the simplest theories for thermal fluctuations in membrane stacks is the "discrete 
harmonic model" (58). It describes membranes without surface tension, and assumes that the fluc- 
tuations are governed by two factors: The bending stiffness K c of single membranes, and a penalty 
for compressing or swelling the stack, which is characterized by a compressibility modulus B. The 
free energy is given in harmonic approximation 

^ DH = £ J/xdy{^(^ + ^f + ^(h n ~h n+1+ df}, (1) 

where d is the average distance between layers. The first part describes the elasticity of single 
membranes, and the second part accounts for the interactions between membranes. Being quadratic, 
the free energy functional (Q~|i is simple enough that statistical averages can be calculated exactly. 

We have tested this theory in detail with coarse-grained molecular simulations (C. Loison et 
al 11591 ) of a binary mixture of amphiphiles and solvent, within a Smit-type model originally intro- 
duced by Soddemann et al |60|. The elementary units are spheres with a hard core radius a, which 
may have two types: "hydrophilic" or "hydrophobic". Beads of the same type attract each other at 
distances less than 1.5a. "Amphiphiles" are tetramers made of two hydrophilic and two hydropho- 
bic beads, and "solvent" particles are single hydrophilic beads. We studied systems of up to 153600 
beads with constant pressure molecular dynamics, using a simulation box of fluctuating size and 
shape, and a Langevin thermostat to maintain constant temperature. 
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Figure 9: Left: Snapshot of a bilayer stack (30720 amphiphiles and 30720 solvent beads). The 
"hydrophobic" tail beads are dark, the "hydrophilic" head beads and the solvent beads are light; 
Right: Snapshot of a single membrane position h n (x, y). 



At suitable pressures and temperatures, the system assumes a fluid lamellar phase. Our configu- 
rations contained up to 15 stacked bilayers, which contained 20 volume percent solvent. A configu- 
ration snapshot is shown in Fig.|9](left). In order to test the theory ([T), one must first determine the 
local positions of the membranes in the stack. This was done as follows: 

1. The simulation box was divided into N x N y N z cells. (N x = N y = 32). The size of the cells 
fluctuates because the dimensions of the box fluctuate. 

2. In each cell, the relative density of tail beads p M (x, y, z) was calculated. It is defined as the 
ratio of the number of tail beads and the total number of beads. 

3. The hydrophobic space was defined as the space where the relative density of tail beads is 
higher than a given threshold po- (The value of the threshold was roughly 0.7, i.e., 80 % of 
the maximum value of p M1 ). 

4. Cells belonging to the hydrophobic space are connected to clusters. Two hydrophobic cells 
that share at least one vortex are attributed to the same cluster. Each cluster defines a mem- 
brane. 

5. For each membrane n and each position (x. y), the two heights ftjf y) and h™"(x, y), where 
the density p m (x, y, z) passes through the threshold po, are determined. The mean position is 
defined as the average h n (x, y) = (h™ n + ft,™")/2. 

The algorithm identifies membranes even if they have pores. At the presence of other defects, such 
as necks or passages (connections between membranes), additional steps must be taken, which are 
not described here. A typical membrane conformation h n (x, y) is shown in Fig. [9] (right). 

The statistical distribution of h n (x,y) can be analyzed and compared with theoretical predic- 
tions. For example, the transmembrane structure factor, which describes correlations between mem- 
brane positions in different membranes is given by 



>(<?) = {h m {q)*h m +n(<l)) = Ml) 



X 



(2) 
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where h n (q) is the Fourier transform of h n (x, y) in the (x, z/)-plane and X = q A K c /B is a dimen- 
sionless parameter. The function So(q) can also be calculated explicitly |59|. The prediction (f2]l 
can be tested by simply plotting the ratio s n /so vs. q for different n. The functional form of the 
curves should be given by the expression in the square brackets, with only one fit parameter K c /B. 
Figure [10] shows the simulation data. The agreement with the theory is very good over the whole 
range of wave vectors q. 

Thus the molecular simulations confirm the validity of the mesoscopic model ([T). Moreover, the 
analysis described above yields a value for the phenomenological parameter K c /B. By analyzing 
other quantities, one can also determine K c and B separately. This gives the elastic parameters of 
membranes and their effective interactions, and establishes a link between the molecular and the 
mesoscopic description. 

Many important characteristics of membranes not only depend on their elasticity, but also on 
their defects. For example, the permeability is crucially influenced by the number and properties 
of membrane pores. A number of atomistic and coarse grained simulation studies have therefore 
addressed pore formation [61 62, 63, 64, 65 1, mostly in membranes under tension. In contrast, the 
membranes in the simulations of Loison et al have almost zero surface tension. This turns out to 
affect the characteristics of the pores quite dramatically Il66ll67l . 

FigureQT]shows a snapshot of a hydrophobic layer which contains a number of pores. The pores 
have nucleated spontaneously. They "live" for a while, grow and shrink without diffusing too much, 
until they finally disappear. Most pores close very quickly, but a few large ones stay open for a long 
time. A closer analysis shows that pores are hydrophilic, i.e., the amphiphiles rearrange themselves 
at the pore edge, such that solvent beads in the pore center are mainly exposed to head beads. The 
total number of pores is rather high in our system. This is because the amphiphiles are short and the 
membranes are thin. 

The analysis algorithm presented above not only localizes membranes, but also identifies pores, 
their position and their shape. Therefore one can again test the appropriate mesoscopic theories for 
pore formation. The simplest Ansatz for the free energy of a pore with the area A and the contour 
length c has the form [68 1 

E = E + Xc--/A, (3) 

where E is a core energy, A a material parameter called line tension, and 7 the surface tension. The 
second term describes the energy penalty at the pore rim. The last term accounts for the reduction 
of energy due to the release of surface tension in a stretched membrane. In our case, the surface 
tension is close to zero (7 sa 0), because the simulation was conducted such that the pressure tensor 
is isotropic. Thus the last term vanishes. 
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Figure 10: Ratios of transmembrane structure factors si/sq and s^j 'sq vs. in-plane wavevector q 
in units of cr _1 . The solid lines correspond to a theoretical fit to Eq. (f2]i with one (common) fit 
parameter K c /B. After Ref. ||59l . 
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Figure 11: Snapshot of a single bilayer (top view. Only hydrophobic beads are shown). 




Figure 12: Distribution of pore contours in a semi-logarithmic plot. Left: Raw data, Right: Divided 
by the degeneracy function g(c). After Ref. [66|. 



In this simple free energy model, the pore shapes should be distributed according to a Boltzmann 
distribution, P(c) oc exp(— Ac). Figure [12] (left) shows a histogram of contour lengths P(c). The 
bare data do not reflect the expected exponential behavior. Something is missing. Indeed, a closer 
look reveals that the naive exponential Ansatz disregards an important effect: The "free energy" 

gives only local free energy contributions, i.e., those stemming from local interactions and local 
amphiphile rearrangements. In addition, one must also account for the global entropy of possible 
contour length conformations. Therefore, we have to evaluate the "degeneracy" of contour lengths 
g(c), and test the relation 

P(c) oc g(c) exp(-Ac). (4) 

Figure [T2] (right) demonstrates that this second Ansatz describes the data very well. From the linear 
fit to the data, one can extract a value for the line tension A. 

The model (01 makes a second important prediction: Since the free energy only depends on the 
contour length, pores with the same contour length are equivalent, and the shapes of these pores 
should be distributed like those of two-dimensional self-avoiding ring polymers. In other words, 
they are not round, but have a fractal structure. From polymer theory, one knows that the size R g of 
a two dimensional self-avoiding polymer scales roughly like R g oc N 3 / 4 with the polymer length 
N, In our case, the "polymer length" is the contour length c. Thus the area A of a pore should scale 
like 

Ao,R 2 g0 ,{C 3 / 4 f ^C 3 ' 2 . (5) 

Figure Qj] shows that this is indeed the case. 

Other properties of pores have been investigated, e.g., dynamical properties, pore lifetimes, pore 
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Figure 13: Pore area vs. contour length (arbitrary units). After Ref. [66 1. 



correlations etc. [66|. All results were in good agreement with the line tension model ©. In sum, 
we find that the fluctuations and (pore) defects in membrane stacks can be described very well 
by a combination of two simple mesoscopic theories: An effective interface model for membrane 
undulations (Eq. (Q]l), and a line tension model for the pores (Eq. (0 |67|. 

The simulations of Loison et al thus demonstrate nicely that simulations of coarse-grained 
molecular models can be used to test mesoscopic theories, and to bridge between the microscopic 
and the mesoscopic level. 

2.3 Mesoscopic membrane models 

When it comes to large scale structures and long time dynamics, simulations of molecular models 
become too cumbersome. It then proves useful to drop the notion of particles altogether and work 
directly with continuous mesoscopic theories. In general, continuum theories are mostly formulated 
in terms of continuous fields, which stand for local, coarse-grained averages of some appropriate 
microscopic quantities. In the case of membranes systems, it is often more adequate to keep the 
membranes as separate objects, despite their microscopic thickness, and treat them as interfaces in 
space. This leads to random-interface theories. 

2.3.1 Random interface models 

In random-interface models, membranes are represented by two-dimensional sheets in space, whose 
conformations are distributed according to an effective interface free energy functional. An example 
of such a functional has already been given in Eq. ([TJ. However, the free energy (Q]l can only be used 
to describe slightly fluctuating planar membranes in the (xy) -plane. In order to treat membranes of 
arbitrary shape (vesicles etc.), one must resort to the formalism of differential geometry [ 12|. 

One of the simplest approaches of this kind is the spontaneous curvature model. The elastic free 
energy of a membrane is given by l70ll 



where JdS- integrates over the surface of the membrane, and H and K are the mean and Gaus- 
sian curvature, which are related to the local curvature radii i?i and i?2 via H = {l/Ri + 1/ i?2)/2, 
K = \/(RiR 2 ) [69]. The parameters k, R (the bending stiffnesses) and H (the spontaneous cur- 
vature) depend on the specific material properties of the membrane. The total membrane surface 
A = JdS is taken to be constant. 

In the case Ho = and k > 0, k > 0, the free energy (O describes membranes that want to be 
flat, and imposes a bending penalty on all local curvatures. H n ^ implies that the membrane has 
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Figure 14: Connectivity update in the dynamic triangulation algorithm 



a favorite mean curvature. It is worth noting that the integral over the Gaussian curvature, JK dS, 
depends solely on the topology of the interfaces. For membranes without rims, the Gauss-Bonnet 
theorem states that 



where the Euler characteristics \e — 2(c — g) counts the number of closed surfaces c (including 
cavities) minus the number of handles g. 

Other free energy expressions have been proposed, which allow for an asymmetric distribution 
of lipids between both sides of the membrane, and/or for fluctuations of A IT7T1I721 . Here, we shall 
only discuss the spontaneous curvature model. Our membranes shall be self-avoiding, i.e., they have 
excluded volume interactions and cannot cross one another. 

Our task is to formulate a simulation model that mimicks the continuous space theory as closely 
as possible. More precisely, we first need a method that samples random self-avoiding surfaces, and 
second a way to discretize the free energy functional (O on those surfaces. 

A widely used method to generate self-avoiding surfaces is the dynamic triangulation algo- 
rithm |[73l[74l [751176117711781 . Surfaces are modeled by triangular networks of spherical beads, which 
are linked by tethers of some given maximum length Iq (tethered-bead model). The tethers define 
the neighbor relations on the surface. In order to enforce self-avoidance, the beads are equipped 
with hard cores (diameter cr), and the maximum tether length is chosen in the range a < Iq < y/3a. 
The surfaces are sampled with Monte Carlo, with two basic types of moves: Moves that change 
the position of beads, and moves that change the connectivity of the network, i.e., redistribute the 
tethers between the chains. The latter is done by randomly cutting tethers and flipping them between 
two possible diagonals of neighboring triangles (Fig. IT4T>. Updates of bead positions are subject to 
the constraint that the maximum tether length must not exceed Iq. In addition, connectivity updates 
have to comply with the requirement that at least three tethers are attached to one bead, and that two 
beads cannot be connected to one another by more than one tether. Otherwise, attempted updates are 
accepted or rejected according to the standard Metropolis prescription. Of course, additional, more 
sophisticated Monte Carlo moves can be introduced as well. 

Next, the free energy functional (|6]l must be adapted for the triangulated surface. If no rims are 
present, we only need an expression for the first term. Since the integral over the Gaussian curvature 
K is then essentially given by the Euler-characteristic, it is much more convenient and accurate to 
evaluate the latter directly (i.e., count closed surfaces minus handles) than to attempt a discretization 



The question how to discretize the bending free energy in an optimal way is highly non-trivial, 
and several approximations have been proposed over the years 1T781 . Here we give a relatively recent 
expression due to Kumar et al in 2001 [79|: The surface integral J dS is replaced by a sum over 
all triangles t with area A t . For each triangle t, one considers the three neighbor triangles t' and 
determines the length L u > of the common side as well as the angle Q t v between the two triangle 




(7) 



of K. 
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Figure 15: Osmotically driven vesicle translocation through a pore. Configuration snapshots. From 
G. T. Linke l87l 

normal vectors. The free energy is then approximated by 

= \ / dA(2H - H f « I £ A t (^^ - tfo) 2 . (8) 

This completes the tethered-bead representation for random interfaces. We have presented one of 
the most simple variants. The models can be extended in various directions, e.g., to represent mixed 
membranes [80|, membranes that undergo liquid/gel transitions [81], membranes with fluctuating 
geometry 11821 or pores J83 1 etc. They have even been used to study vesicle dynamics, e.g., vesicles in 
shear flow [84, 85] . The dynamics can be made more realistic by using hybrid algorithms, where the 
beads are moved according to molecular dynamics or Brownian dynamics, and the tethers are flipped 
with Monte Carlo. Unfortunately, the tether flips cannot easily be integrated in a pure molecular 
dynamics simulation. This is a slight drawback of the tethered-bead approach. 

An interesting alternative has recently been proposed by Noguchi et al [86]. He developed 
a specific type of potential, which forces (single) beads to self-assemble into a two-dimensional 
membrane without being connected by tethers. The membrane is coarse-grained as a curved surface 
as in tethered-bead models, but tethers are no longer necessary. This opens the way for a fully 
consistent treatment of dynamics in membrane systems with frequent topology changes. 

2.3.2 Application example: passage of vesicles through pores 

As an example for an application of a tethered-bead model, we shall discuss a problem from the 
physics of vesicles. Linke, Lipowsky and Gruhn have recently investigated the question whether a 
vesicle can be forced to cross a narrow pore by osmotic pressure gradients l87l[88l . Vesicle passage 
through pores plays a key role in many strategies for drug delivery. On passing through the pore, the 
vesicle undergoes a conformational change, which is expensive and creates a barrier. On the other 
hand, the vesicle can reduce its energy, if it crosses into a region with lower osmotic pressure. Linke 
et al studied the question, whether the barrier height can be reduced by a sufficient amount that the 
pore is crossed spontaneously. 

To this end, they considered a vesicle filled with N osmotically active particles. The concentra- 
tion of these particles outside the vesicle is given by C12 on the two sides of the pore. For a vesicle in 
the process of crossing the pore from side 1 to side 2, the vesicle area on both sides of the membrane 
is denoted Ai 2 , and the volumes are V\ ,2- Hence the osmotic contribution to the free energy is 

T„ = (-JVln((7i + V 2 )/V ) + ciVi + C2V2), (9) 

where Vq is some reference volume. In the Monte Carlo simulations, a series of different crossing 
stages was sampled (see Fig. IT5J. This was done by introducing an appropriate reweighting factor, 
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which kept the area A2 within a certain, pre-determined range. Several simulations were done with 
different windows for A%. The procedure has similarity with umbrella sampling, except that the 
windows did not overlap. Nevertheless, the extrapolation of the histograms allowed to evaluate the 
free energy quite accurately, as a function of A2 . 

The resulting free energy landscape is shown in Fig. [16] for different osmolarities c\ at the side 
where the vesicle starts. As expected, one finds a free energy barrier for low c\. On increasing c\, 
the height of the barrier decreases, until it finally disappears. Hence the vesicle can be forced to cross 
the pore spontaneously, if c\ is chosen sufficiently high. In contrast, changing the bending rigidity 
has hardly any effect on the height of the potential barrier. 
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Figure 16: Free energy landscape on crossing the pore for different osmolarities c\ starting from 
c\ = 2700 (open circles) to c\ — 3100 (plus) in steps of 100. All units are arbitrary. The osmolarity 
C2 is kept fixed. See text for more explanation. FromG. Linkeefa/|87 88 1. 

This application shows how mesoscopic models of membrane systems can be used to study the 
physical basis of a potentially technologically relevant process. A comparable study with a particle- 
based model would have been much more expensive, and both the setup and the analysis would 
have been more difficult: One has no well-defined interface, one has only indirect information on 
the bending stiffness etc. The simplifications in the model clearly help not only to get results more 
quickly, but also to understand them more systematically. 

2.4 Summary 

In this section, we have introduced and discussed two important classes of coarse-grained models 
for membrane systems: Bead-spring models for simulations on the coarse-grained molecular level, 
and random interface models for the coarse-grained mesoscopic level. We have shown that such 
studies can give insight into generic properties of membranes and into basic physical mechanisms 
on different scales. We have also shown that it is possible to bridge between different levels, and 
that simulations of particle-based models can be used to test mesoscopic theories. 

3 Liquid crystals and surfactant layers under shear 

So far, we have talked about equilibrium simulations. The special properties of complex fluids 
lead to particularly peculiar behavior at nonequilibrium. In this section, we will discuss coarse- 
grained simulations of complex fluids under shear. As in the previous section, we will move from 
the coarse-grained molecular level to the mesoscopic level. However, we will now focus on the 
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simulation methods rather than on the construction of simulation models. As example systems, we 
will consider surfactant layers and liquid crystals. 

3.1 Introduction 

3.1.1 Strain rate and shear stress 

We begin with recalling some basic fluid mechanics ||89ll . In the continuum description, the fluid is 
thought to be divided into fluid elements, which are big enough that microscopic details are washed 
out, but small enough that each element can still be considered to be homogoneous. On a macro- 
scopic scale, the fluid is described by spatially varying fields, e.g., the velocity or flow field u(r, t). 
The velocity gradient in the direction perpendicular to the flow is called the local strain rate. 

One of the central quantities in fluid mechanics is the stress tensor a. It describes the forces 
that the surrounding fluid exerts on the surfaces of a fluid element: For a surface with orientation 
N (\N\ — 1), the force per area is Fi/A = J^j a ij^j- The diagonal components of a give the 
longitudinal stress and correspond to normal forces. The off-diagonal components give the shear 
stress and correspond to tangential forces. 

A fluid element with the velocity u is thus accelerated according to 

i 

where p is the local mass density. The stress tensor dij is often divided in three parts. 

(Tij = -pS ij +a^ ie + a' ij (11) 

The first term describes the effect of the local pressure. It is always present, even in an equilibrium 
fluid at rest. The second term, off**, accounts for elastic restoring forces, if applicable. The last 
term, •, gives the contribution of viscous forces and is commonly called the viscous stress tensor. 
It vanishes in the absence of velocity gradients, i.e., if diVj = 0. The concrete relation between the 
stress tensor and the local fields is called constitutive equation and characterizes the specific fluid 
material. 

In standard hydrodynamics, one considers so-called Newtonian fluids, which have no elasticity 
^eiasiic = q-) an( j a ]j near5 instantaneous relation between the viscous stress tensor and the velocity 
gradient, 

cr'ij = ^ Li J kl dkUl - ( 12 ) 

kl 

In the case of an incompressible, isotropic, Newtonian fluid, one has only one independent material 
parameter, the kinematic viscosity rj, and the constitutive equation reads 

(Tij = -P Sij + rj(d t u :j + djUi), (13) 

where the value of the local pressure p(r) is determined by the requirement that the incompressibility 
condition 

W=0. (14) 

is always fulfilled. Equations (ITOb in combination with ( fT3l are the famous Navier-stokes equations 
(without external forces). 

Complex fluids have internal degrees of freedom, and large characteristic length and time scales. 
Therefore, they become non-Newtonian already at low shear stress, and standard hydrodynamics 
does not apply. The relation between the stress tensor and the velocity gradient is commonly not 
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Figure 17: Sketches of nonlinear stress-strain flow curves Left: Shear thickening and shear thinning; 
Right: Nonmonotonic stress-strain relation leading to phase separation. The horizontal dotted lines 
connect the two coexisting phases for three selected global strain rates. 

linear, sometimes one has no unique relation at all. If the characteristic time scales are large, one 
encounters memory effects, e.g., the viscosity may decrease with time under shear stress l|2l |4l [89l . 

These phenomena are often considered in planar Couette geometry: The flow u points in the 
x-direction, a velocity gradient in the y-direction is imposed, and the relation between the strain rate 
7 = du x /dy and the shear stress a xy is investigated. For Newtonian fluids, the two are proportional. 
In complex fluids, the differential viscosity da xy /d^i may decrease with increasing strain rate (shear 
thinning), or increase (shear thickening). In some materials, the stress-strain flow curve may even be 
nonmonotonic. Since such flow curves are mechanically unstable in macroscopic systems, this leads 
to phase separation, and the system becomes inhomogeneous (Fig. \Yl\ [ 90 9l1 l92ll93ll94ll95ll96l 
I97ll98ll99l 11031 [T04l . The separation of a fluid into two states with distinctly different strain rates 
is called shear banding. As in equilibrium phase transitions, the signature of the coexistence region 
is a plateau region in the resulting stress-strain flow curve of the inhomogeneous systems. However, 
one has two important differences to the equilibrium case: First, the location of the plateau cannot 
be determined by a Maxwell construction, but has to be calculated by solving the full dynamical 
equations. Second, the plateau is not necessarily horizontal. Even though mechanical stability 
requires that two coexisting phases must have the same shear stress, the "plateau" may have a slope 
or even be curved because different coexisting phases may be selected for different shear rates [100 
[1051 (see Fig. [171 

3.1.2 Nematic liquid crystals 

Next we recapitulate some basic facts on liquid crystals. The building blocks of these materials are 
anisotropic particles, e.g., elongated molecules, associated structures such as wormlike micelles, or 
even rodlike virusses. A liquid crystal phase is an intermediate state between isotropic liquid and 
crystalline solid, which is anisotropic (particles are orientationally ordered), but in some sense still 
fluid (particles are spatially disordered in at least one direction). Some such structures are sketched 
in Fig. [18] The phase transitions are in some cases triggered by the temperature (thermotropic liquid 
crystals), and in other cases by the concentration of the anisotropic particles (lyotropic liquid crys- 
tals). Lamellar membrane stacks (see previous section) are examples of smectic liquid crystalline 
structures. In this section, we shall mainly be concerned with nematics, where the particles are 
aligned along one common direction, but otherwise fluid. 

For symmetry reasons, the transition between the isotropic fluid and the nematic fluid must be 
first order. The nematic order is commonly characterized by the symmetric traceless order tensor 



Qij = (-^fididj - 1)), 



(15) 
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where d is a unit vector characterizing the orientation of a particle. In an isotropic fluid, Q is 
zero. Otherwise, the largest eigenvalue gives the nematic order parameter S, and the corresponding 
eigenvector is the direction of preferred alignment, the "director" n. Note that only the direction of 
n matters, not the orientation, i.e., the vectors n and — n are equivalent. 

Since the ordered nematic state is degenerate with respect to a continuous symmetry (the di- 
rection of the director), it exhibits elasticity: The free energy costs of spatial director variations 
must vanish if the wavelength of the modulations tends to infinity. At finite wavelength, the sys- 
tem responds to such distortions with elastic restoring forces. Symmetry arguments show that these 
depend only on three independent material constants 11161 , the Frank constants Ku, K22, and ^33. 
The elastic free energy is given by [ 14 1 

J 7 — = JdfiKuiVrt) 2 + K 22 (n • (V x n)f + K 33 (n x (V x n)) 2 }. (16) 

Figure [T9l illustrates the corresponding fundamental distortions, the splay {Ku), twist (K22), and 
bend mode (X33). 

How does shear affect a liquid crystalline fluid? In the isotropic phase and at low strain rate , the 
situation is still comparably simple. The velocity gradient imposes a background angular velocity 

n=I(Vxn) (17) 

on the fluid, which the particles pick up 11061 . Since they rotate faster while perpendicular to the 
flow, they also acquire a very slight effective orientation at the angle it/ A relative to the flow direc- 
tion. At higher strain rate, the order parameter increases, and the alignment angle decreases 11071 . 
This sometimes even leads to a shear-induced transition into the nematic phase 110811109111 lOlfTTTl 

Deep in the nematic phase, the nematic fluid can be described by the flow field u(r, t) and the 
director field n(r,t) (\n\ = 1). We shall now assume incompressibility (VS = 0) and a linear, 
instantaneous relation stress-strain relation as in Eq. ( TT2b . Compared to simple fluids, there are 
several complications. 

First, nematic fluids are elastic, hence the stress tensor <7y ( fTTT l has an elastic contribution 

^V— '>!"■■ ■ (18) 

J 8{d l n k ) 

The free energy !F amic is given by Eq. dT6t . 
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Figure 18: Examples of liquid crystalline phases. In the nematic phase (N), the particles are oriented, 
but translationally disordered. In the smectic phases, they form layers in one direction, but remain 
fluid within the layers. In the smectic B Hai and smectic F phase, the structure within the layers is not 
entirely disordered, but has a type of order called hexatic. See Refs. 1 14 15] for explanation. The 
isotropic phase (I) is entirely disordered (regular isotropic liquid). 
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Second, the fluid is locally anisotropic, therefore the constitutive equation for the viscous stress 
tensor contains more independent terms than in simple fluids. The structure of a'^ which is compat- 
ible with all symmetry requirements reads |fT31l 



a ij = ai X! n i n j n kniA k i + a 2 riiNj + a 3 NiUj (19) 

kl 

+ ai + a 5 ^2 turikA-jk + ae ^ njUkAik, 

k k 

where A is the symmetric flow deformation tensor, Aij — \{diUj + djUi), and the vector N gives 
the rate of change of the director relative to the angular velocity of the fluid, N — 4| — Cl x fl. The 
parameters ai are the so-called Leslie coefficients. They are linked by one relation |fT31 . 

a2 + a-3 = ae — as (20) 

hence one has five independent material constants. 

Third, the flow field and the director field are coupled. Thus one must take into account the 
dynamical evolution of the director field, 

pihj = gj + diTTij, (21) 

which is driven by the intrinsic body force g, and the director stress tensor Wij. Here p\ is the 
moment of inertia per unit volume. As for the stress tensor, one can derive the general form of the 
constitutive equations for g and 7r with symmetry arguments. The individual expressions can be 
found in the literature ifTSIl . Here, we only give the final total dynamical equation for n, 

Pl n 3 = An, + £ d z (^g-y ) - g - 7l iV- ■ - 72 £ mAij , (22) 

where 71 = a 3 — a 2 , 72 = «6 — a 5, and A is a Lagrange parameter which keeps |n| constant. 

To summarize, the dynamical equations of nematohydrodynamics for incompressible nematic 
fluids are given by Eqs. dTob and (1221 together with (fTTT i. dT6l >. ( fT8l ). and dT9b . These are the Leslie- 
Ericksen equations of nemato-hydrodynamics Q3) . The theory depends on five independent vis- 
cosity coefficients ai and on three independent elastic constants Ku. Not surprisingly, it predicts a 
rich spectrum of phenomena. For example, in the absence of director gradients, steady shear flow is 
found to have two opposing effects on the director: It rotates the molecules, and it aligns them. The 
relative strength of these two tendencies depends on the material parameters. As a result, one may 
either encounter stable flow alignment or a rotating "tumbling" state II 15111 16lll 171 [TT8l . 

The Ericksen-Leslie theory describes incompressible and fully ordered nematic fluids. The basic 
version presented here does not account for the possibility of disclination lines and other topological 
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Figure 19: Elastic modes in nematic liquid crystals. 
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defects, and it does not include the coupling with a variable density and/or order parameter field. 
The situation is even more complicated in the vicinity of stable or metastable (hidden) thermody- 
namic phase transitions. Thermodynamic forces may contribute to a nonmonotonic stress-strain 
relationship, which eventually lead to mechanical phase separation |98 1. 

Figure [20] shows an experimental example of a nonequilibrium phase diagram for a system of 
rodlike virusses 11021 . 
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Figure 20: Nonequilibrium phase diagram of a lyotropic liquid crystal under shear as a function 
of strain rate 7 and particle density C. (C ncm is the density of the equilibrium isotropic/nematic 
transition). The system is a mixture of rodlike virusses (fd virus) and polymers (dextran). From 

Ref. una. 



3.2 Simulating shear on the particle level: NEMD 

In the linear response regime, i.e., in the low-shear limit, many rheological properties of materials 
can be determined from equilibrium simulations. For example, the viscosity coefficients introduced 
in the previous section can be related to equilibrium fluctuations with Green-Kubo relations Ml 191 

nsa. 

However, the non-Newtonian character of most complex fluids and their resulting unique prop- 
erties manifest themselves only beyond the linear response regime. To study these, non-equilibrium 
simulations are necessary. We will now briefly review non-equilibrium molecular dynamics (NEMD) 
methods for the simulation of molecular model systems under shear. Some of the methods are cov- 
ered in much more detail in Refs. II 1 2 1 1 1 1 22l 1 1231 . 

In NEMD simulations, one faces two challenges. First, one must mechanically impose shear. 
Second, most methods that enforce shear constantly pump energy into the system. Hence one must 
get rid of that heat, i.e., apply an appropriate thermostat. We will address these two issues separately. 

The most direct way of imposing shear is to confine the system between two rough walls, and 
either move one of them (for Couette flow), or apply a pressure gradient (for Poiseuille flow). Varnik 
and Binder have shown that this approach can be used to measure the shear viscosity in polymer 
melts [124|. It has the merit of being physical: Strain is enforced physically, and the heat can 
be removed in a physical way by coupling a thermostat to the walls. On the other hand, the system 
contains two surfaces, and depending on the material under consideration, one may encounter strong 
surface effects. This can cause problems. 
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An alternative straightforward way to generate planar Couette flow is to use moving periodic 
boundary conditions as illustrated in Fig. [21] (Lees-Edwards boundary conditions 01211 11251 ): In 
order to enforce shear flow u — (u x ,0, 0) with an average strain rate 7 = du x /dy, one proceeds 
as follows: One replicates the particles in the x and z direction like in regular periodic boundary 
conditions, 



v — ► v 



x.y 7 ' x.y 

r z -> r z ±L z . (23) 



In the y-direction, the replicated particles acquire an additional velocity U = jL y and an offset U t 
in the x-direction, 

r x — > r x ± (Ut modL x ) 
r y ~ T r y i L y 
r z -> r z 
v — > v±Ue x 

Here L x , ytZ are the dimensions of the simulation box, and e x the unit vector in the a;-direction. 

Lees-Edwards boundary conditions are perfectly sufficient to enforce planar Couette flow. In 
some applications, it has nevertheless proven useful to supplement them with fictitious bulk forces 
that favor a linear flow profile. One particularly popular algorithm of this kind is the SLLOD algo- 
rithm [122, 1261. For the imposed flow u(r) — jye x , the SLLOD equations of motion for particles 
i read 



dn _ p 
dt 

%_ 

dt 



lVie x (25) 

mi 

= Fi-jp'^, (26) 



where = rrii(vi — u(rl)) is the momentum of particle i in a reference frame moving with the local 
flow velocity u(ri), and Fj is the regular force acting on the particle i. The equations of motion 
can be integrated with standard techniques, or with specially devised operator-split algorithms 111281 



11291 . It is also possible to formulate SLLOD variants for arbitrary steady flow 111271 . The SLLOD 
algorithm is very useful for the determination of viscosities in complex fluids 11301 11 191 11311 11321 
11331 11341 . It has the advantage that the system relaxes faster towards a steady state, and that the 
analysis of the data is simplified I122L 

A fourth method to enforce shear has recently been proposed by Miiller-Plathe [ 135|. One uses 
regular periodic boundary conditions, and divides the system into slabs in the desired direction of 
flow gradient. In periodic intervals, the particle in the central slab (y = 0) with the largest velocity 
component in the x direction, and the particle in the topmost slab (y — L y /2) with the largest 




Figure 21: Lees-Edwards boundary conditions. 
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Figure 22: Schematic sketch of the Muller-Plathe method for imposing shear. In periodical intervals, 
the particle in the middle slab that moves fastest in the x-direction and the particle in the top slab that 
moves fastest in the (— a;)-direction exchange their momentum components p x . After Ref. [ 135 1. 



velocity component in the (—x) -direction is determined, and the momentum components p x of the 
two particles are swapped. This leads to a zigzag-shaped shear profile (Fig.l22l. 

The Muller-Plathe algorithm can be implemented very easily. Just one small change turns a 
molecular dynamics program for thermal equilibrium into a NEMD program that produces shear 
flow. Furthermore, the algorithm conserves momentum and energy: Contrary to the algorithms 
discussed before, it does not pump heat into the system. Instead, the heat is constantly redistributed. 
It is drained out of the system at y — and y = L y /2 by the momentum swaps, but it is also 
produced in the bulk through the energy dissipation associated with the shear flow. As a result, the 
zigzag shear profile is accompanied with a W-shaped temperature profile. Since this is not always 
wanted, the algorithm is sometimes used in combination with a thermostat. A slight drawback of the 
algorithm is the fact that it produces intrinsically inhomogeneous profiles, due to the presence of the 
two special slabs at y = and y = L y /2. This may cause problems in small systems. 

Having introduced these four methods to impose shear, we now turn to the problem of ther- 
mostatting. All shearing methods except for the Muller-Plathe method produce heat, therefore the 
thermostat is an essential part of the simulation algorithm. 

A thermostat is a piece of algorithm that manipulates the particle velocities such that the tem- 
perature is adjusted to its desired value. This is commonly done by adjusting the kinetic energy. 
In a flowing fluid with flow velocity u(f), the kinetic energy has a flow contribution and a thermal 
contribution. Therefore, it is important to define the thermostat in a frame that moves with the fluid, 
i.e., to couple it to the "peculiar" velocities = Vi — u{fi), rather than to the absolute velocities 
Vi. Unfortunately, most thermostats don't account for this automatically, and one has to put in the 
flow profile by hand. The only exception is the dissipative particle dynamics (DPD) thermostat (see 
below). 

This raises the question how to determine the flow u{r). In a homogeneous system, at low strain 
rates, it can be acceptable to use the pre-imposed profile (u(r) — jye x in our earlier examples). 
Thermostats that use the pre-imposed flow profile are called "profile-biased thermostats". At high 
strain rates, and in inhomogeneous systems, the use of profile-biased thermostats is dangerous and 
may twist the results 1122111331 . In that case, the profiles u(r) must be calculated directly from the 
simulations, "on the fly". Such thermostats are called "profile-unbiased thermostats". 

One can distinguish between deterministic and stochastic thermostatting methods. The simplest 
deterministic thermostat simply rescales the peculiar velocities of all particles after every time step 
such that the total thermal kinetic energy remains constant. Among the more sophisticated determin- 
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istic thermostats, the ones that have been used most widely in NEMD simulations are the so-called 
Gaussian thermostats 0122L and the Nose Hoover thermostats 021111361 . 

The idea of Gaussian thermostats is very simple: For a system of particles j subject to forces 

fj, one minimizes the quantity C = J2j(Pj ~ fj) 2 /(^ m j) w i m respect to pp with the constraint 
that either the time derivative of the total energy, or the time derivative of the total kinetic energy 
be zero (Gaussian isoenergetic and Gaussian isothermal thermostat, respectively). In the absence 
of constraints, the minimization procedure would yield Newton's law. The constraints modify the 
equations of motion in a way that can be considered in some sense minimal. One obtains 

Pj=fj-otfj, (27) 

where a is a Lagrange multiplier, which has to be chosen such that the desired constraint - constant 
energy or constant kinetic energy - is fulfilled. 

The Nose Hoover thermostats are well-known from equilibrium molecular dynamics and shall 
not be discussed in detail here. We only note that in NEMD simulations, one uses not only the 
standard isothermal Nose-Hoover thermostats, but also isoenergetic variants. 

Deterministic thermostats have the slightly disturbing feature that they change the dynamical 
equations in a rather unphysical way. It is sometimes hard to say how that might affect the system. 
At equilibrium, Nose-Hoover thermostats rest on a well-established theoretical basis: They can be 
derived from extended Hamiltonians, they produce the correct thermal averages for static quantities 
etc. Outside of equilibrium, the situation is much more vague. For the case of steady states, some 
results on the equivalence of thermostats are fortunately available. Evans and Sarman have shown 
that steady state averages and time correlation functions are identical for Gaussian isothermal and 
isoenergetic thermostats and for Nose-Hoover thermostats 111231 11371 . Ruelle has recently proved 
that the Gaussian isothermal and the Gaussian isoenergetic thermostat are equivalent in an infinite 
system which is ergodic under space translations 0381 . 

An alternative approach to thermostatting are the stochastic methods, such as the Langevin ther- 
mostat II 13911 and the DPD thermostat II 1401 . They maintain the desired temperature by introducing 
friction and stochastic forces. This approach has the virtue of being physically motivated. The dissi- 
pative and stochastic terms stand for microscopic degrees of freedom, which have been integrated out 
in the coarse-grained model, but can still absorb and carry heat. From the theory of coarse-graining, 
it is well-known that integrating out degrees of freedom effectively turns deterministic equations of 
motion into stochastic equations of motion 11 1 4 11 11421 11431 . 

The Langevin thermostat is the simplest stochastic dynamical model. One modifies the equation 
of motion for the peculiar momenta by introducing a dissipative friction £ and a random force ffj, 

(28) 

where ffj fulfills 

(rfk) = (29) 
(WO^CO) = 2(k B TS l3 5 a(3 6(t-t'), (30) 

i.e., the noise is delta-correlated and Gaussian distributed with mean zero and variance 2(ksT 
(fluctuation-dissipation theorem). In practice, this can be done by adding {~C v 'j a ) an< ^ a random 
number to each force component fi a in each time step At. The random number must be distributed 
with mean zero and variance y/2ftBT(/At. Because of the central limit theorem, the distribution 
does not necessarily have to be Gaussian, if the time step is sufficiently small 1 144|. 

The DPD thermostat is an application of the dissipative particle dynamics algorithm 111401 . The 
idea is similar to that of the Langevin algorithm, but instead of coupling to absolute velocities, the 
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thermostat couples to velocity differences between neighbor particles. Therefore, the algorithm is 
Galilean invariant, and accounts automatically for the difference between flow velocity and thermal 
velocity. We refer to B. Dunweg's lecture (this book) for an introduction into the DPD method (see 
also E6]). 

This completes our brief introduction into NEMD methods for systems under shear. We close 
with a general comment: The strain rate 7 has the dimension of inverse time. It introduces a time 
scale I/7 in the system, which diverges in the limit 7—^0. Therefore, systems at low strain rates 
converge very slowly towards their final steady state, and the study of systems at low strain rates is 
very time-consuming. In fact, even in coarse-grained molecular models, simulated strain rates are 
usually much higher than experimentally accessible strain rates. 

We now give two examples of recent large-scale NEMD studies of inhomogeneous complex 
systems under shear. 

3.2.1 First application example: Nematic-Isotropic Interfaces under shear 

Our first example is a study of the behavior of a Nematic-Isotropic interface under shear (Germano 
and Schmid II 1451 11461 ). Interfaces play a central role for shear banding, and understanding their 
structure should provide a key to the understanding of phase separation under shear. Our study is 
a first step in this direction. We will discuss the questions whether the structure of an equilibrium 
interface is affected by shear, whether the phase transition is shifted, and whether shear banding can 
be observed. 

The model system was a fluid of soft repulsive ellipsoids with the aspect ratio 15:1, in a simula- 
tion box of size (150:300:150) (in units of the ellipsoid diameter a). The density was chosen in the 
coexistence region between the nematic and isotropic phase, such that the system phase separates 
at equilibrium. The initial configuration was a relaxed equilibrium configuration 2.2with a nematic 
slab and an isotropic slab, separated by two interfaces. The latter align the director of the nematic 
phase such that it is parallel to the interface l|1471ll48l[T49l . 




Figure 23: Snapshot of a sheared nematic-isotropic interface (1 15200 particles) at strain rate 7 = 
O.OOIt -1 . Also shown is the coordinate system used throughout this section. From Ref. II 1451 - 
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Figure 24: Illustration of block analysis to obtain local profiles. 



Shear was imposed with Lees-Edwards boundary conditions in combination with a profile un- 
biased Nose Hoover thermostat, in the direction normal to the interface. The two coexisting phases 
thus share the same shear stress. Two setups are possible within this "common stress" geometry. 
In the flow-aligned setup, the director points in the direction of flow u (the x direction); in the log- 
rolling setup, it points in the direction of "vorticity" V x u (the z direction). Only the flow-aligned 
setup turned out to be stable in our system. We considered mainly the strain rate 7 = O.OOlr -1 . 
Here r is the natural time unit, t = a^Jm/kBT, with the particle mass m, the particle diame- 
ter a, and the temperature T. This strain rate is small enough that the interface is still stable. A 
configuration snapshot is shown in Fig. [23] 

To analyze the system, it was split into columns of size B x L y x B. The columns were further 
split into bins, which contained enough particles that a local order parameter could be determined. In 
this manner, one obtains local order parameter profiles for each column in each configuration. From 
the profile S(y), we determined the local positions of the two interfaces as follows: We computed at 
least two coarse-grained profiles S(y; Si) by convoluting the profile S(y) with a symmetrical box- 
like coarse-graining function of width 5,-. The coarse-graining width must be chosen such that the 
coarse-grained profiles still reflect the interfaces, but short ranged fluctuations are averaged out. The 
intersection of two averaged profiles locates a "dividing surface", where the negative order parameter 
excess on the nematic side balances the positive order parameter excess on the paranematic side. We 
define this to be the "local position" of the interface. Due to fluctuations, it depends slightly on the 
particular choice of the <5;. Nevertheless, the procedure gives remarkably unambiguous values even 
for strongly fluctuating profiles II 14911 . The procedure is illustrated in Fig.l24l 

Once the positions tiNi,hiN of the two interfaces have been determined, we calculate profiles 
for all quantities of interest and shift them by the amount h^i or Kin, respectively. This allows to 
perform averages over local profiles. Furthermore, the interface positions hjN and h^i themselves 
can be used to analyze the fluctuation spectrum of the interface positions. 
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Figure 25: Interfacial order parameter profile at strain rate 7 = O.OOIt -1 (flow aligned setup). Inset 
shows the difference between profiles with and without shear. From Ref. 1 146 1. 
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Unexpectedly, the shear turned out to have almost no influence at all on the interfacial structure. 
Under shear, the interfaces are slightly broadened, indicating that the interfacial tension might be 
slightly reduced (Fig. \25[ . But the effect is barely noticeable. Likewise, the density profile and 
the capillary wave spectrum are almost undistinguishable from those of the equilibrium interface, 
within the error. The velocity profiles, however, are much more interesting. The flow profile exhibits 
a distinct kink at the interface position (Fig. [26]). Hence we clearly observe shear banding - the total 
strain is distributed nonuniformly between the two phases. At the same shear stress, the nematic 
phase supports a higher strain rate than the isotropic phase. Second and surprisingly, the interface 
apparently also induces a streaming velocity gradient in the vorticity direction (the z-direction). As 
a consequence, vorticity flow is generated in opposite directions in the isotropic and in the nematic 
phase. This flow is symmetry breaking, and as yet, we have no explanation for it. We hope that 
future theoretical and simulation studies will clarify this phenomenon. 

At high strain rates, the interface gets destroyed. Figure [27] shows a nonequilibrium phase dia- 
gram, which has been obtained from interface simulations of small systems. The local strain rate is 
always higher in the nematic phase than in the isotropic phase. It is remarkable that the coexistence 
region does not close up. Instead, the interface disappears abruptly at average strain rates above 
7 w 0.006/t. 
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Figure 27: Phase diagram of the flow-aligned system. The dashed tie lines connect coexisting states. 
I denotes the isotropic region, and N the nematic region. From Ref. 1 146 1. 



3.2.2 Second application example: Shear-induced phenomena in surfactant layers 

As a second example, we discuss the effect of shear on lamellar stacks of surfactants. This has been 
studied in coarse-grained molecular simulations by Soddemann, Guo, Kremer et al II150II1511 . The 
model system was very similar to that used in Section l2.2.4l except that it does not contain solvent 
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Figure 28: Reorientation of surfactant layers under shear from an initial transverse state at low strain 
rate, 7 = 0.001/r. Reprinted with permission from Ref. |150|. Copyright 2002 by the American 
Physical Society. 
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Figure 29: Reorientation of surfactant layers under shear from an initial transverse state at high strain 
rate, 7 = 0.03/t. Reprinted with permission from Ref. |150|. Copyright 2002 by the American 
Physical Society. 



particles ll60l . Shear was induced with the Muller-Plathe algorithm II 1 351 in combination with a 
DPP thermostat l60l. 

The first series of simulations by Guo et al 015011 addressed the question, whether the surfactant 
layers can be forced to reorient under shear. Three different orientations of the layers with respect 
to the shear geometry were considered: In the transverse orientation, the layers are normal to the 
flow direction, in the parallel orientation, they are normal to the direction of the velocity gradient, 
and in the perpendicular orientation, they lie in the plane of the flow and the velocity gradient. 
The transverse orientation is unstable under shear flow. Experimentally, both a transition from the 
transverse to the parallel orientation and from the transverse to the perpendicular orientation have 
been observed lTT32l[T53l . 

Guo et al have studied this by simulations of a system of dimer amphiphiles. Figures|28land|29l 
show a series of snapshots for different times at two different strain rates. Both systems were set 
up in the transverse state. At low strain rate, the shear generates defects in the transverse lamellae. 
Mediated by these defects, the lamellae gradually reorient into the parallel state. At high strain rate, 
the lamellae are first completely destroyed and then reorganize in the perpendicular state. 

The final state is not necessarily the most favorable steady state. In fact, Guo et al found that 
the strain energy dissipation was always smallest in the perpendicular state, regardless of the strain 
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Figure 30: Shear induced undulations. Left: configuration snapshot at strain rate 7 = 0.015/t. 
Right: Undulation amplitude as a function of strain rate. Undulations set in at a strain rate of 
7 fa 0.01. From Ref. fT5fl . 
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Figure 3 1 : Snapshots from a Lattice-Boltzmann simulation of a sheared hybrid aligned nematic cell. 
White regions are ordered, and black regions are disordered. The numbers give time in milliseconds. 
FromRef. lfT60l . 

rate. When starting from the transverse state at low strain rates, the parallel state formed nevertheless 
for kinetic reasons. At intermediate strain rates, shear may induce undulations in the parallel state. 
This has been predicted by Auernhammer et al Ml 541 11551 and studied by Soddemann et al by 
computer simulation of a system of layered tetramers H 15 lH . The phenomenon results from the 
coupling between the layer normal, the tilt of the molecules, and the shear flow. It is triggered by 
the fact that shear flow induces tilt. The tilted surfactant layer dilates and uses more area, which 
eventually leads to an undulation instability. This could indeed be observed in the simulations. 
Figure[30]shows a snapshot of an undulated configuration, and a plot of the undulation amplitude as 
a function of strain rate. The undulations set in at a well-defined strain-rate. 

Our two examples show that coarse-grained molecular simulations of complex, inhomogeneous 
fluids under shear are now becoming possible. A wealth of new, intriguing physics can be expected 
from this field in the future. 

3.3 Simulations at the Mesoscopic Level 

We will only touch very briefly on the wide field of mesoscopic simulations for complex fluids under 
shear. 

The challenge of mesoscopic simulations is to find and/or formulate the appropriate mesoscopic 
model, and then put it on the computer. For liquid crystals, we have already discussed one candidate, 
the Leslie-Ericksen theory. However, this is by far not the only available mesoscopic theory for liquid 
crystals. A popular alternative for the description of lyotropic liquid crystals is the Doi theory, which 
is based on a Smoluchowski equation for the distribution function for rods. Numerous variants and 
extensions of the Doi model have been proposed. A relatively recent review on both the Leslie- 
Ericksen and the Doi theory can be found in Ref . H156II . 

The phenomenological equations can be solved with traditional methods of computational fluid 
dynamics. An alternative approach which is particularly well suited to simulate flows in complex and 
fluctuating geometries is the Lattice Boltzmann method, which has been discussed by B. Diinweg 
in his lecture. Yeomans and coworkers have recently proposed a Lattice Boltzmann algorithm for 
liquid-crystal hydrodynamics, which allows to simulate liquid crystal flow 11571 [1581. The start- 
ing point is the Beris-Edwards theory 115911 . a dynamical model for liquid crystals that consists of 
coupled dynamical equations for the velocity profiles and the order tensor Q (Eq. fTBY The Lattice- 
Boltzmann scheme works with the usual discrete velocities i and distributions fi(r n ) for the fluid 
flow field on the lattice site f n . In addition, a second distribution Gj(r„) is introduced, which de- 
scribes the flow of the order tensor field. The time evolution includes free streaming and collision 
steps. The exact equations are quite complicated and shall not be given here, they can be found in 
Refs. l[T58l . 

As an application example of this method, we discuss a recent simulation by Marenduzzo et al 
of a sheared hybrid aligned nematic (HAN) cell 0160111611 . The surfaces in a HAN cell impose con- 
flicting orientations on the director of an adjacent nematic fluid, i.e., one surface orients ("anchors") 
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the fluid in a parallel way ("planar"), and the other in a perpendicular way ("homeotropic"). In the 
system under consideration, the parameters were chosen such that the fluid is in the isotropic phase, 
but very close to the nematic phase. The fluid is at rest at time t > 0. At times t > 0, the surfaces 
are moved relative to each other, and drag the fluid with them (no-slip boundary conditions). Figure 
[5T| shows snapshots of the cell for different times. The evolution of the structure inside the cell is 
quite complex. First, two ordered bands form close to the surface, due to the fact that the flow field 
still builds up and strain rates are quite high in the vicinity of the surfaces. The director in these 
bands is flow-aligned, and in the lower band, its direction competes with the anchoring energy of the 
homeotropic surface. As a result, the lower band unbinds and crosses the system to join the other 
band. Finally, the top band also unbinds and moves towards the center of the cell. This simulation 
shows nicely how a complex dynamical process, which results from an interplay of shear flow and 
elastic deformations, can be studied with a mesoscopic simulation method. We note that the time 
scale of this simulation is seconds, i.e., inaccessible for molecular simulations. 

3.4 Summary 

Under shear, complex fluids exhibit a wealth of new, interesting, and practically relevant phenom- 
ena. To study these, various simulation approaches for different levels of coarse-graining have been 
developed. In this section, we have presented and discussed several variants of non-equilibrium 
molecular dynamics (NEMD), and illustrated their use with examples of large scale NEMD simula- 
tions of inhomogeneous liquid crystals and surfactant layers. Furthermore, we have briefly touched 
on mesoscopic simulation methods for liquid crystals under shear. 

4 Conclusions 

Due to the rapidly improving performance of modern computers, complex fluids can be studied 
on a much higher level today than just ten years ago. A large number of coarse-grained models 
and methods have been designed, which allow to investigate different aspects of these materials, 
at equilibrium as well as far from equilibrium. The new computational possibilities have vitally 
contributed to the current boost of soft matter science. 

In the two central sections of this lecture, we have given an introduction into two important 
chapters of computational soft matter: In Section |2] we have given a rough overview over idealized 
coarse-grained simulation models for membranes, and hopefully conveyed an idea of the poten- 
tial of such models. In Section [3] we have discussed simulation methods for complex fluids under 
shear. Generally, nonequilibrium studies of complex fluids are still much more scarce than equilib- 
rium studies. Much remains to be done, both from the point of view of method development and 
applications, and many exciting developments can be expected for the near future. 
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